/*
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE aomatrix_test
#include <boost/test/unit_test.hpp>
#include <votca/xtp/aomatrix.h>
#include <votca/xtp/orbitals.h>

using namespace votca::xtp;
using namespace votca;
using namespace std;

BOOST_AUTO_TEST_SUITE(aomatrix_test)

QMMolecule Methane() {
  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();
  QMMolecule mol(" ", 0);
  mol.LoadFromFile("molecule.xyz");
  return mol;
}

BOOST_AUTO_TEST_CASE(aomatrices_test) {

  ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << endl;
  basisfile << "  <element name=\"H\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "  <element name=\"C\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "</basis>" << endl;
  basisfile.close();

  QMMolecule mol = Methane();
  BasisSet basis;
  basis.Load("3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, mol);
  AOOverlap overlap;
  overlap.Fill(aobasis);
  Eigen::MatrixXd overlap_ref = Eigen::MatrixXd::Zero(17, 17);

  overlap_ref << 1, 0.191448, 0, 0, 0, 0.180314, 0, 0, 0, 0.0189724, 0.0808612,
      0.0189724, 0.0808612, 0.0189724, 0.0808612, 0.0189724, 0.0808612,
      0.191448, 1.00001, 0, 0, 0, 0.761361, 0, 0, 0, 0.194748, 0.401447,
      0.194748, 0.401447, 0.194748, 0.401447, 0.194748, 0.401447, 0, 0, 1, 0, 0,
      0, 0.528959, 0, 0, 0.169584, 0.135615, 0.169584, 0.135615, -0.169584,
      -0.135615, -0.169584, -0.135615, 0, 0, 0, 1, 0, 0, 0, 0.528959, 0,
      0.169584, 0.135615, -0.169584, -0.135615, -0.169584, -0.135615, 0.169584,
      0.135615, 0, 0, 0, 0, 1, 0, 0, 0, 0.528959, 0.169584, 0.135615, -0.169584,
      -0.135615, 0.169584, 0.135615, -0.169584, -0.135615, 0.180314, 0.761361,
      0, 0, 0, 1, 0, 0, 0, 0.338796, 0.668849, 0.338796, 0.668849, 0.338796,
      0.668849, 0.338796, 0.668849, 0, 0, 0.528959, 0, 0, 0, 1, 0, 0, 0.290649,
      0.340149, 0.290649, 0.340149, -0.290649, -0.340149, -0.290649, -0.340149,
      0, 0, 0, 0.528959, 0, 0, 0, 1, 0, 0.290649, 0.340149, -0.290649,
      -0.340149, -0.290649, -0.340149, 0.290649, 0.340149, 0, 0, 0, 0, 0.528959,
      0, 0, 0, 1, 0.290649, 0.340149, -0.290649, -0.340149, 0.290649, 0.340149,
      -0.290649, -0.340149, 0.0189724, 0.194748, 0.169584, 0.169584, 0.169584,
      0.338796, 0.290649, 0.290649, 0.290649, 1, 0.645899, 0.00778321, 0.116994,
      0.00778321, 0.116994, 0.00778321, 0.116994, 0.0808612, 0.401447, 0.135615,
      0.135615, 0.135615, 0.668849, 0.340149, 0.340149, 0.340149, 0.645899, 1,
      0.116994, 0.354983, 0.116994, 0.354983, 0.116994, 0.354983, 0.0189724,
      0.194748, 0.169584, -0.169584, -0.169584, 0.338796, 0.290649, -0.290649,
      -0.290649, 0.00778321, 0.116994, 1, 0.645899, 0.00778321, 0.116994,
      0.00778321, 0.116994, 0.0808612, 0.401447, 0.135615, -0.135615, -0.135615,
      0.668849, 0.340149, -0.340149, -0.340149, 0.116994, 0.354983, 0.645899, 1,
      0.116994, 0.354983, 0.116994, 0.354983, 0.0189724, 0.194748, -0.169584,
      -0.169584, 0.169584, 0.338796, -0.290649, -0.290649, 0.290649, 0.00778321,
      0.116994, 0.00778321, 0.116994, 1, 0.645899, 0.00778321, 0.116994,
      0.0808612, 0.401447, -0.135615, -0.135615, 0.135615, 0.668849, -0.340149,
      -0.340149, 0.340149, 0.116994, 0.354983, 0.116994, 0.354983, 0.645899, 1,
      0.116994, 0.354983, 0.0189724, 0.194748, -0.169584, 0.169584, -0.169584,
      0.338796, -0.290649, 0.290649, -0.290649, 0.00778321, 0.116994,
      0.00778321, 0.116994, 0.00778321, 0.116994, 1, 0.645899, 0.0808612,
      0.401447, -0.135615, 0.135615, -0.135615, 0.668849, -0.340149, 0.340149,
      -0.340149, 0.116994, 0.354983, 0.116994, 0.354983, 0.116994, 0.354983,
      0.645899, 1;

  bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.0001);
  BOOST_CHECK_EQUAL(check_overlap, 1);
  if (!check_overlap) {
    cout << "ref" << endl;
    cout << overlap_ref << endl;
    cout << "result" << endl;
    cout << overlap.Matrix() << endl;
  }

  AOKinetic kinetic;
  kinetic.Fill(aobasis);
  Eigen::MatrixXd kinetic_ref = Eigen::MatrixXd::Zero(17, 17);
  kinetic_ref << 16.579, -1.43503, 0, 0, 0, 0.10275, 0, 0, 0, -0.0439437,
      0.0214514, -0.0439437, 0.0214514, -0.0439437, 0.0214514, -0.0439437,
      0.0214514, -1.43503, 1.35738, 0, 0, 0, 0.346414, 0, 0, 0, -0.0154143,
      0.103305, -0.0154143, 0.103305, -0.0154143, 0.103305, -0.0154143,
      0.103305, 0, 0, 2.58667, 0, 0, 0, 0.41751, 0, 0, 0.0928809, 0.0755731,
      0.0928809, 0.0755731, -0.0928809, -0.0755731, -0.0928809, -0.0755731, 0,
      0, 0, 2.58667, 0, 0, 0, 0.41751, 0, 0.0928809, 0.0755731, -0.0928809,
      -0.0755731, -0.0928809, -0.0755731, 0.0928809, 0.0755731, 0, 0, 0, 0,
      2.58667, 0, 0, 0, 0.41751, 0.0928809, 0.0755731, -0.0928809, -0.0755731,
      0.0928809, 0.0755731, -0.0928809, -0.0755731, 0.10275, 0.346414, 0, 0, 0,
      0.293786, 0, 0, 0, 0.0889197, 0.139112, 0.0889197, 0.139112, 0.0889197,
      0.139112, 0.0889197, 0.139112, 0, 0, 0.41751, 0, 0, 0, 0.489642, 0, 0,
      0.169257, 0.135141, 0.169257, 0.135141, -0.169257, -0.135141, -0.169257,
      -0.135141, 0, 0, 0, 0.41751, 0, 0, 0, 0.489642, 0, 0.169257, 0.135141,
      -0.169257, -0.135141, -0.169257, -0.135141, 0.169257, 0.135141, 0, 0, 0,
      0, 0.41751, 0, 0, 0, 0.489642, 0.169257, 0.135141, -0.169257, -0.135141,
      0.169257, 0.135141, -0.169257, -0.135141, -0.0439437, -0.0154143,
      0.0928809, 0.0928809, 0.0928809, 0.0889197, 0.169257, 0.169257, 0.169257,
      1.5494, 0.293152, -0.0206172, -0.00736852, -0.0206172, -0.00736852,
      -0.0206172, -0.00736852, 0.0214514, 0.103305, 0.0755731, 0.0755731,
      0.0755731, 0.139112, 0.135141, 0.135141, 0.135141, 0.293152, 0.274788,
      -0.00736852, 0.0301943, -0.00736852, 0.0301943, -0.00736852, 0.0301943,
      -0.0439437, -0.0154143, 0.0928809, -0.0928809, -0.0928809, 0.0889197,
      0.169257, -0.169257, -0.169257, -0.0206172, -0.00736852, 1.5494, 0.293152,
      -0.0206172, -0.00736852, -0.0206172, -0.00736852, 0.0214514, 0.103305,
      0.0755731, -0.0755731, -0.0755731, 0.139112, 0.135141, -0.135141,
      -0.135141, -0.00736852, 0.0301943, 0.293152, 0.274788, -0.00736852,
      0.0301943, -0.00736852, 0.0301943, -0.0439437, -0.0154143, -0.0928809,
      -0.0928809, 0.0928809, 0.0889197, -0.169257, -0.169257, 0.169257,
      -0.0206172, -0.00736852, -0.0206172, -0.00736852, 1.5494, 0.293152,
      -0.0206172, -0.00736852, 0.0214514, 0.103305, -0.0755731, -0.0755731,
      0.0755731, 0.139112, -0.135141, -0.135141, 0.135141, -0.00736852,
      0.0301943, -0.00736852, 0.0301943, 0.293152, 0.274788, -0.00736852,
      0.0301943, -0.0439437, -0.0154143, -0.0928809, 0.0928809, -0.0928809,
      0.0889197, -0.169257, 0.169257, -0.169257, -0.0206172, -0.00736852,
      -0.0206172, -0.00736852, -0.0206172, -0.00736852, 1.5494, 0.293152,
      0.0214514, 0.103305, -0.0755731, 0.0755731, -0.0755731, 0.139112,
      -0.135141, 0.135141, -0.135141, -0.00736852, 0.0301943, -0.00736852,
      0.0301943, -0.00736852, 0.0301943, 0.293152, 0.274788;

  bool check_kinetic = kinetic.Matrix().isApprox(kinetic_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_kinetic, 1);
  if (!check_kinetic) {
    cout << "ref" << endl;
    cout << kinetic_ref << endl;
    cout << "result" << endl;
    cout << kinetic.Matrix() << endl;
  }

  AOCoulomb coulomb;
  coulomb.Fill(aobasis);
  Eigen::MatrixXd coulomb_ref = Eigen::MatrixXd::Zero(17, 17);
  coulomb_ref << 1.66592, 4.0152, 0, 0, 0, 5.96515, 0, 0, 0, 1.86584, 4.83598,
      1.86584, 4.83598, 1.86584, 4.83598, 1.86584, 4.83598, 4.0152, 18.379, 0,
      0, 0, 31.3743, 0, 0, 0, 10.3031, 26.6091, 10.3031, 26.6091, 10.3031,
      26.6091, 10.3031, 26.6091, 0, 0, 4.75942, 0, 0, 0, 7.02939, 0, 0, 1.96588,
      2.48722, 1.96588, 2.48722, -1.96588, -2.48722, -1.96588, -2.48722, 0, 0,
      0, 4.75942, 0, 0, 0, 7.02939, 0, 1.96588, 2.48722, -1.96588, -2.48722,
      -1.96588, -2.48722, 1.96588, 2.48722, 0, 0, 0, 0, 4.75942, 0, 0, 0,
      7.02939, 1.96588, 2.48722, -1.96588, -2.48722, 1.96588, 2.48722, -1.96588,
      -2.48722, 5.96515, 31.3743, 0, 0, 0, 64.1609, 0, 0, 0, 21.3671, 58.4239,
      21.3671, 58.4239, 21.3671, 58.4239, 21.3671, 58.4239, 0, 0, 7.02939, 0, 0,
      0, 21.387, 0, 0, 5.07919, 8.88641, 5.07919, 8.88641, -5.07919, -8.88641,
      -5.07919, -8.88641, 0, 0, 0, 7.02939, 0, 0, 0, 21.387, 0, 5.07919,
      8.88641, -5.07919, -8.88641, -5.07919, -8.88641, 5.07919, 8.88641, 0, 0,
      0, 0, 7.02939, 0, 0, 0, 21.387, 5.07919, 8.88641, -5.07919, -8.88641,
      5.07919, 8.88641, -5.07919, -8.88641, 1.86584, 10.3031, 1.96588, 1.96588,
      1.96588, 21.3671, 5.07919, 5.07919, 5.07919, 13.9085, 26.8615, 5.54669,
      17.0411, 5.54669, 17.0411, 5.54669, 17.0411, 4.83598, 26.6091, 2.48722,
      2.48722, 2.48722, 58.4239, 8.88641, 8.88641, 8.88641, 26.8615, 68.5967,
      17.0411, 50.7702, 17.0411, 50.7702, 17.0411, 50.7702, 1.86584, 10.3031,
      1.96588, -1.96588, -1.96588, 21.3671, 5.07919, -5.07919, -5.07919,
      5.54669, 17.0411, 13.9085, 26.8615, 5.54669, 17.0411, 5.54669, 17.0411,
      4.83598, 26.6091, 2.48722, -2.48722, -2.48722, 58.4239, 8.88641, -8.88641,
      -8.88641, 17.0411, 50.7702, 26.8615, 68.5967, 17.0411, 50.7702, 17.0411,
      50.7702, 1.86584, 10.3031, -1.96588, -1.96588, 1.96588, 21.3671, -5.07919,
      -5.07919, 5.07919, 5.54669, 17.0411, 5.54669, 17.0411, 13.9085, 26.8615,
      5.54669, 17.0411, 4.83598, 26.6091, -2.48722, -2.48722, 2.48722, 58.4239,
      -8.88641, -8.88641, 8.88641, 17.0411, 50.7702, 17.0411, 50.7702, 26.8615,
      68.5967, 17.0411, 50.7702, 1.86584, 10.3031, -1.96588, 1.96588, -1.96588,
      21.3671, -5.07919, 5.07919, -5.07919, 5.54669, 17.0411, 5.54669, 17.0411,
      5.54669, 17.0411, 13.9085, 26.8615, 4.83598, 26.6091, -2.48722, 2.48722,
      -2.48722, 58.4239, -8.88641, 8.88641, -8.88641, 17.0411, 50.7702, 17.0411,
      50.7702, 17.0411, 50.7702, 26.8615, 68.5967;
  bool check_coulomb = coulomb.Matrix().isApprox(coulomb_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_coulomb, 1);
  if (!check_coulomb) {
    cout << "ref" << endl;
    cout << coulomb_ref << endl;
    cout << "result" << endl;
    cout << coulomb.Matrix() << endl;
  }

  Eigen::MatrixXd ps_invSqrt = coulomb.Pseudo_InvSqrt(1e-7);
  Eigen::MatrixXd coulombinvsqrt_ref = Eigen::MatrixXd::Zero(17, 17);
  coulombinvsqrt_ref << 1.16163, -0.307153, 9.97835e-16, -1.10548e-16,
      1.00508e-16, 0.0492417, -1.93564e-15, 1.85948e-16, -1.33024e-16,
      0.0226533, -0.00856968, 0.0226533, -0.00856968, 0.0226533, -0.00856968,
      0.0226533, -0.00856968, -0.307153, 0.829138, -1.96961e-15, 3.70613e-16,
      -3.72916e-16, -0.608979, 3.60427e-15, -5.74697e-16, 5.22424e-16,
      -0.0218795, 0.0808469, -0.0218795, 0.0808469, -0.0218795, 0.0808469,
      -0.0218795, 0.0808469, 9.97835e-16, -1.96961e-15, 0.679881, 2.75893e-16,
      3.26907e-17, 2.17372e-15, -0.207947, -3.5126e-16, -4.68835e-16,
      -0.0646361, 0.0554385, -0.0646361, 0.0554385, 0.0646361, -0.0554385,
      0.0646361, -0.0554385, -1.10548e-16, 3.70613e-16, 2.76025e-16, 0.679881,
      -9.91215e-17, -6.96261e-16, 8.15636e-17, -0.207947, 3.66591e-16,
      -0.0646361, 0.0554385, 0.0646361, -0.0554385, 0.0646361, -0.0554385,
      -0.0646361, 0.0554385, 1.00508e-16, -3.72916e-16, 3.28582e-17,
      -9.70689e-17, 0.679881, 5.26717e-16, 1.97761e-16, 5.74318e-16, -0.207947,
      -0.0646361, 0.0554385, 0.0646361, -0.0554385, -0.0646361, 0.0554385,
      0.0646361, -0.0554385, 0.0492417, -0.608979, 2.17372e-15, -6.96261e-16,
      5.26717e-16, 1.36693, -3.70734e-15, 9.17285e-16, -8.05922e-16, -0.119143,
      -0.233467, -0.119143, -0.233467, -0.119143, -0.233467, -0.119143,
      -0.233467, -1.93564e-15, 3.60427e-15, -0.207947, 8.14087e-17, 1.96659e-16,
      -3.70734e-15, 0.526054, 1.22817e-16, 2.41716e-16, -0.012717, -0.162726,
      -0.012717, -0.162726, 0.012717, 0.162726, 0.012717, 0.162726, 1.85948e-16,
      -5.74697e-16, -3.50753e-16, -0.207947, 5.70762e-16, 9.17285e-16,
      1.1573e-16, 0.526054, -7.70072e-16, -0.012717, -0.162726, 0.012717,
      0.162726, 0.012717, 0.162726, -0.012717, -0.162726, -1.33024e-16,
      5.22424e-16, -4.6736e-16, 3.63935e-16, -0.207947, -8.05922e-16,
      2.42516e-16, -7.66731e-16, 0.526054, -0.012717, -0.162726, 0.012717,
      0.162726, -0.012717, -0.162726, 0.012717, 0.162726, 0.0226533, -0.0218795,
      -0.0646361, -0.0646361, -0.0646361, -0.119143, -0.012717, -0.012717,
      -0.012717, 0.606742, -0.164568, 0.0157438, 0.0286141, 0.0157438,
      0.0286141, 0.0157438, 0.0286141, -0.00856968, 0.0808469, 0.0554385,
      0.0554385, 0.0554385, -0.233467, -0.162726, -0.162726, -0.162726,
      -0.164568, 0.512794, 0.0286141, -0.0727664, 0.0286141, -0.0727664,
      0.0286141, -0.0727664, 0.0226533, -0.0218795, -0.0646361, 0.0646361,
      0.0646361, -0.119143, -0.012717, 0.012717, 0.012717, 0.0157438, 0.0286141,
      0.606742, -0.164568, 0.0157438, 0.0286141, 0.0157438, 0.0286141,
      -0.00856968, 0.0808469, 0.0554385, -0.0554385, -0.0554385, -0.233467,
      -0.162726, 0.162726, 0.162726, 0.0286141, -0.0727664, -0.164568, 0.512794,
      0.0286141, -0.0727664, 0.0286141, -0.0727664, 0.0226533, -0.0218795,
      0.0646361, 0.0646361, -0.0646361, -0.119143, 0.012717, 0.012717,
      -0.012717, 0.0157438, 0.0286141, 0.0157438, 0.0286141, 0.606742,
      -0.164568, 0.0157438, 0.0286141, -0.00856968, 0.0808469, -0.0554385,
      -0.0554385, 0.0554385, -0.233467, 0.162726, 0.162726, -0.162726,
      0.0286141, -0.0727664, 0.0286141, -0.0727664, -0.164568, 0.512794,
      0.0286141, -0.0727664, 0.0226533, -0.0218795, 0.0646361, -0.0646361,
      0.0646361, -0.119143, 0.012717, -0.012717, 0.012717, 0.0157438, 0.0286141,
      0.0157438, 0.0286141, 0.0157438, 0.0286141, 0.606742, -0.164568,
      -0.00856968, 0.0808469, -0.0554385, 0.0554385, -0.0554385, -0.233467,
      0.162726, -0.162726, 0.162726, 0.0286141, -0.0727664, 0.0286141,
      -0.0727664, 0.0286141, -0.0727664, -0.164568, 0.512794;
  bool check_coulombinvsqrt = ps_invSqrt.isApprox(coulombinvsqrt_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_coulombinvsqrt, 1);
  if (!check_coulombinvsqrt) {
    cout << "ref" << endl;
    cout << coulombinvsqrt_ref << endl;
    cout << "result" << endl;
    cout << ps_invSqrt << endl;
  }

  Eigen::MatrixXd ps_invSqrtgw = coulomb.Pseudo_InvSqrt_GWBSE(overlap, 1e-7);
  Eigen::MatrixXd coulombinvsqrtgw_ref = Eigen::MatrixXd::Zero(17, 17);
  coulombinvsqrtgw_ref << 1.18216, -0.221181, 3.31007e-16, -2.77017e-17,
      5.4648e-17, -0.024762, -4.67451e-17, 2.42791e-18, -7.39467e-17, 0.0191022,
      0.00100653, 0.0191022, 0.00100653, 0.0191022, 0.00100653, 0.0191022,
      0.00100653, -0.387336, 0.905821, -9.70665e-16, 1.12468e-16, -7.20284e-17,
      -0.458094, -2.51973e-16, 2.48351e-16, 4.71999e-16, 0.00887867, 0.00238888,
      0.00887867, 0.00238888, 0.00887867, 0.00238888, 0.00887867, 0.00238888,
      4.23087e-16, -6.09524e-16, 0.705798, 4.81316e-16, 1.3367e-16, 3.47315e-17,
      -0.13188, -1.27084e-16, 3.48185e-17, -0.0613791, 0.0311124, -0.0613791,
      0.0311124, 0.0613791, -0.0311124, 0.0613791, -0.0311124, 1.0619e-17,
      4.2506e-17, 4.78159e-16, 0.705798, -1.81479e-16, 7.98158e-17,
      -7.45323e-17, -0.13188, 1.77653e-16, -0.0613791, 0.0311124, 0.0613791,
      -0.0311124, 0.0613791, -0.0311124, -0.0613791, 0.0311124, 8.8482e-17,
      -1.10927e-16, 7.25119e-17, -1.27572e-16, 0.705798, 9.94042e-17,
      1.74247e-16, 9.21421e-17, -0.13188, -0.0613791, 0.0311124, 0.0613791,
      -0.0311124, -0.0613791, 0.0311124, 0.0613791, -0.0311124, 0.140772,
      -0.845289, 1.08307e-15, -1.38689e-16, 2.19957e-17, 1.29026, 6.23914e-16,
      -5.71606e-16, -7.72672e-16, -0.162867, -0.0533785, -0.162867, -0.0533785,
      -0.162867, -0.0533785, -0.162867, -0.0533785, -2.36379e-16, 1.27991e-16,
      -0.281562, -1.88063e-16, 1.68519e-16, 4.23859e-16, 0.519419, 2.15099e-16,
      -2.20189e-16, -0.00395714, -0.139109, -0.00395714, -0.139109, 0.00395714,
      0.139109, 0.00395714, 0.139109, -2.23714e-16, 4.63617e-16, -3.5789e-16,
      -0.281562, 1.691e-16, -6.9076e-16, 4.32912e-16, 0.519419, -1.50764e-16,
      -0.00395714, -0.139109, 0.00395714, 0.139109, 0.00395714, 0.139109,
      -0.00395714, -0.139109, -1.29103e-16, 2.23031e-16, -1.79135e-16,
      1.09785e-16, -0.281562, -3.18918e-16, -2.12058e-16, 3.09506e-17, 0.519419,
      -0.00395714, -0.139109, 0.00395714, 0.139109, -0.00395714, -0.139109,
      0.00395714, 0.139109, 0.0210071, -0.0124102, -0.0704645, -0.0704645,
      -0.0704645, -0.0878611, -0.00143061, -0.00143061, -0.00143061, 0.628099,
      -0.0864583, 0.0165705, 0.0141384, 0.0165705, 0.0141384, 0.0165705,
      0.0141384, -0.0234403, 0.126287, 0.0856204, 0.0856204, 0.0856204,
      -0.249203, -0.182259, -0.182259, -0.182259, -0.210269, 0.423453,
      0.0481652, -0.0956113, 0.0481652, -0.0956113, 0.0481652, -0.0956113,
      0.0210071, -0.0124102, -0.0704645, 0.0704645, 0.0704645, -0.0878611,
      -0.00143061, 0.00143061, 0.00143061, 0.0165705, 0.0141384, 0.628099,
      -0.0864583, 0.0165705, 0.0141384, 0.0165705, 0.0141384, -0.0234403,
      0.126287, 0.0856204, -0.0856204, -0.0856204, -0.249203, -0.182259,
      0.182259, 0.182259, 0.0481652, -0.0956113, -0.210269, 0.423453, 0.0481652,
      -0.0956113, 0.0481652, -0.0956113, 0.0210071, -0.0124102, 0.0704645,
      0.0704645, -0.0704645, -0.0878611, 0.00143061, 0.00143061, -0.00143061,
      0.0165705, 0.0141384, 0.0165705, 0.0141384, 0.628099, -0.0864583,
      0.0165705, 0.0141384, -0.0234403, 0.126287, -0.0856204, -0.0856204,
      0.0856204, -0.249203, 0.182259, 0.182259, -0.182259, 0.0481652,
      -0.0956113, 0.0481652, -0.0956113, -0.210269, 0.423453, 0.0481652,
      -0.0956113, 0.0210071, -0.0124102, 0.0704645, -0.0704645, 0.0704645,
      -0.0878611, 0.00143061, -0.00143061, 0.00143061, 0.0165705, 0.0141384,
      0.0165705, 0.0141384, 0.0165705, 0.0141384, 0.628099, -0.0864583,
      -0.0234403, 0.126287, -0.0856204, 0.0856204, -0.0856204, -0.249203,
      0.182259, -0.182259, 0.182259, 0.0481652, -0.0956113, 0.0481652,
      -0.0956113, 0.0481652, -0.0956113, -0.210269, 0.423453;
  bool check_coulombinvsqrtgw =
      ps_invSqrtgw.isApprox(coulombinvsqrtgw_ref, 0.00001);
  BOOST_CHECK_EQUAL(check_coulombinvsqrtgw, 1);
  if (!check_coulombinvsqrtgw) {
    cout << "ref" << endl;
    cout << coulombinvsqrtgw_ref << endl;
    cout << "result" << endl;
    cout << ps_invSqrtgw << endl;
  }
}

BOOST_AUTO_TEST_CASE(aomatrices_contracted_test) {

  std::ofstream basisfile("contracted.xml");
  basisfile << "<basis name=\"cc-pVTZ\">" << std::endl;
  basisfile << "<element name=\"C\">" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "			<constant decay=\"8236.0\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.0005424302\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"1235.0\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.0041964279\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"280.8\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.0215409141\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"79.27\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.0836149496\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"25.59\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.2398716189\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"8.997\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.4437518201\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"3.319\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.3535796965\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"0.3643\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.0091763661\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "			<constant decay=\"8236.0\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.0001963922\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"1235.0\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.0015259503\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"280.8\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.007890449\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"79.27\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.0315148705\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"25.59\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.0969100083\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"8.997\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.2205415263\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"3.319\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"-0.2960691129\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"0.3643\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"1.0405034329\" type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "			<constant decay=\"0.9059\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"S\">" << std::endl;
  basisfile << "			<constant decay=\"0.1285\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"S\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"P\">" << std::endl;
  basisfile << "			<constant decay=\"18.71\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.0394263872\" type=\"P\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"4.133\">"
            << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.2440889849\" type=\"P\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "			<constant decay=\"1.2\">" << std::endl;
  basisfile << "				<contractions "
               "factor=\"0.8154920089\" type=\"P\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"P\">" << std::endl;
  basisfile << "			<constant decay=\"0.3827\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"P\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"P\">" << std::endl;
  basisfile << "			<constant decay=\"0.1209\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"P\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"D\">" << std::endl;
  basisfile << "			<constant decay=\"1.097\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"D\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"D\">" << std::endl;
  basisfile << "			<constant decay=\"0.318\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"D\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "		<shell scale=\"1.0\" type=\"F\">" << std::endl;
  basisfile << "			<constant decay=\"0.761\">"
            << std::endl;
  basisfile << "				<contractions factor=\"1.0\" "
               "type=\"F\"/>"
            << std::endl;
  basisfile << "			</constant>" << std::endl;
  basisfile << "		</shell>" << std::endl;
  basisfile << "	</element>" << std::endl;
  basisfile << "</basis>" << std::endl;
  basisfile.close();

  std::ofstream xyzfile("C.xyz");
  xyzfile << " 1" << std::endl;
  xyzfile << " C" << std::endl;
  xyzfile << " C            .000000     .000000     .000000" << std::endl;
  xyzfile.close();

  QMMolecule mol("C", 0);
  mol.LoadFromFile("C.xyz");
  BasisSet basis;
  basis.Load("contracted.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, mol);
  AOOverlap overlap;
  overlap.Fill(aobasis);
  Eigen::MatrixXd overlap_ref =
      Eigen::MatrixXd::Zero(aobasis.AOBasisSize(), aobasis.AOBasisSize());
  overlap_ref << 1, -0.271918, 0.510892, 0.140478, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.271918, 1, 0.553776,
      0.755963, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0.510892, 0.553776, 1, 0.535796, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.140478, 0.755963,
      0.535796, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.611573, 0, 0, 0.221927, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, -2.7742e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      1, 0, 0, 0.611573, 0, 0, 0.221927, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      2.95039e-17, 0, 0, 0, 5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.611573,
      0, 0, 0.221927, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      -5.55112e-17, 0, 0, 0, 0, 0.611573, 0, 0, 1, 0, 0, 0.674476, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 1.11022e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0.611573, 0, 0, 1, 0, 0, 0.674476, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      -5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.611573, 0, 0, 1, 0, 0,
      0.674476, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.11022e-16, 0, 0, 0, 0, 0,
      0, 0, 0, 0.221927, 0, 0, 0.674476, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 1.11022e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.221927, 0, 0,
      0.674476, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4.16334e-17, 0, 0,
      0, -5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0.221927, 0, 0, 0.674476, 0, 0, 1,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -5.55112e-17, 0, 0, 0, 5.55112e-17, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
      0, 0.531577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.531577, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
      0, 0, 0, 0.531577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0.531577, 0, 0, 0, 0, 1, 0, 0, 0, -5.55112e-17, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0, 1, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0.531577, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.531577,
      -5.55112e-17, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.7742e-17, 0,
      0, 1.11022e-16, 0, 0, 1.11022e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
      0, 0, 0, -5.93751e-18, 0, 0, 0, 0, 0, 0, 0, 2.95039e-17, 0, 0,
      -5.55112e-17, 0, 0, -4.16334e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
      0, 0, 0, 8.37926e-18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.11022e-16, 0, 0,
      -5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 8.37926e-18,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
      0, 0, 0, -5.93751e-18, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5.55112e-17, 0, 0,
      0, 0, 0, -5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.37926e-18, 0,
      0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -5.55112e-17, 0, 0, 0, 0, 0, 5.55112e-17, 0,
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.37926e-18, 0, 0, 0, 1;
  bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.0001);
  if (!check_overlap) {
    std::cout << std::endl;
    std::cout << "Ref" << std::endl;
    std::cout << overlap_ref << std::endl;
    std::cout << "Result" << std::endl;
    std::cout << overlap.Matrix() << std::endl;
  }
  BOOST_CHECK_EQUAL(check_overlap, 1);
}

BOOST_AUTO_TEST_CASE(aocoulomb_inv_test) {
  QMMolecule mol = Methane();
  BasisSet basis;
  basis.Load("3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, mol);

  AOCoulomb cou;
  cou.Fill(aobasis);

  Eigen::MatrixXd PseudoInvSqrt = cou.Pseudo_InvSqrt(1e-7);

  Eigen::MatrixXd Reformed = PseudoInvSqrt * PseudoInvSqrt * cou.Matrix();

  bool check_inv = Reformed.isApprox(Eigen::MatrixXd::Identity(17, 17), 0.0001);
  if (!check_inv) {
    std::cout << "reformed" << endl;
    std::cout << Reformed << endl;
  }
  BOOST_CHECK_EQUAL(check_inv, 1);
}

BOOST_AUTO_TEST_CASE(large_l_test) {
  std::ofstream xyzfile("C2.xyz");
  xyzfile << " 2" << std::endl;
  xyzfile << " C2" << std::endl;
  xyzfile << " C            .000000     .000000     .000000" << std::endl;
  xyzfile << " C            1.000000     .000000     .000000" << std::endl;
  xyzfile.close();

  QMMolecule mol("C", 0);
  mol.LoadFromFile("C2.xyz");

  ofstream basisfile("G.xml");
  basisfile << "<basis name=\"G\">" << endl;
  basisfile << "  <element name=\"C\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"G\">" << endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"G\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "</basis>" << std::endl;
  basisfile.close();

  ofstream auxbasisfile("I.xml");
  auxbasisfile << "<basis name=\"I\">" << endl;
  auxbasisfile << "  <element name=\"C\">" << endl;
  auxbasisfile << "    <shell scale=\"1.0\" type=\"I\">" << endl;
  auxbasisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  auxbasisfile << "        <contractions factor=\"1.562850e-01\" type=\"I\"/>"
               << endl;
  auxbasisfile << "      </constant>" << endl;
  auxbasisfile << "    </shell>" << endl;
  auxbasisfile << "  </element>" << endl;
  auxbasisfile << "</basis>" << std::endl;
  auxbasisfile.close();

  BasisSet basisset;
  basisset.Load("G.xml");

  BasisSet auxbasisset;
  auxbasisset.Load("I.xml");
  AOBasis dftbasis;
  dftbasis.Fill(basisset, mol);

  AOBasis auxbasis;
  auxbasis.Fill(auxbasisset, mol);

  AOOverlap overlap;
  overlap.Fill(auxbasis);

  Index auxbasissize = 26;
  Index dftbasissize = 18;
  Eigen::MatrixXd overlap_ref =
      Eigen::MatrixXd::Zero(auxbasissize, auxbasissize);
  overlap_ref << 1, 0, 0, 0, 0.00980207, 0, 0, 0, 4.12241e-16, 0, 0, 0,
      1.249e-16, 0.0237511, 0, 0, 0, -0.0176617, 0, 0, 0, 0.0148613, 0, 0, 0,
      0.00194066, 0, 1, 0, 0, 0, 2.88549e-16, 0, 0, 0, 3.43734e-16, 0, 0, 0, 0,
      0.00415193, 0, 0, 0, -0.00971203, 0, 0, 0, 0.00833903, 0, 0, 0, 0, 0, 1,
      0, 0, 0, 9.21638e-18, 0, 0, 0, 4.41139e-16, 0, 0, 0, 0, -0.0433016, 0, 0,
      0, 0.0310482, 0, 0, 0, -0.0111976, 0, 0, 0, 0, 0, 1, 0.00473485, 0, 0,
      -8.01766e-17, 0, 0, 0, -2.61143e-16, 0, 0, 0, 0, -0.0138466, -0.00256163,
      0, 0, 0.0185842, 0, 0, 0, -0.00613319, 0, 0.00980207, 0, 0, 0.00473485,
      1.19531, 0, 0, -0.0103735, 0.0155603, 0, 0, 0.0210687, -0.0210687,
      -0.0176617, 0, 0, -0.00256163, 0.05203, 0, 0, 0.00561224, -0.00671312, 0,
      0, -0.0113985, -0.0213299, 0, 2.88549e-16, 0, 0, 0, 1, 0, 0, 0,
      1.2767e-16, 0, 0, 0, 0, -0.00971203, 0, 0, 0, 0.0246417, 0, 0, 0,
      -0.0263528, 0, 0, 0, 0, 0, 9.21638e-18, 0, 0, 0, 1, 0, 0, 0, 4.28107e-17,
      0, 0, 0, 0, 0.0310482, 0, 0, 0, -0.0296094, 0, 0, 0, 0.0227903, 0, 0, 0,
      0, 0, -8.01766e-17, -0.0103735, 0, 0, 1, 0, 0, 0, -6.61878e-17, 0, 0, 0,
      0, 0.0185842, 0.00561224, 0, 0, -0.030451, 0, 0, 0, 0.0211324, 0,
      4.12241e-16, 0, 0, 0, 0.0155603, 0, 0, 0, 1, 0, 0, 0, 1.11022e-16,
      0.0148613, 0, 0, 0, -0.00671312, 0, 0, 0, 0.0281051, 0, 0, 0, -0.0185897,
      0, 3.43734e-16, 0, 0, 0, 1.2767e-16, 0, 0, 0, 1, 0, 0, 0, 0, 0.00833903,
      0, 0, 0, -0.0263528, 0, 0, 0, 0.0451492, 0, 0, 0, 0, 0, 4.41139e-16, 0, 0,
      0, 4.28107e-17, 0, 0, 0, 1, 0, 0, 0, 0, -0.0111976, 0, 0, 0, 0.0227903, 0,
      0, 0, -0.0414569, 0, 0, 0, 0, 0, -2.61143e-16, 0.0210687, 0, 0,
      -6.61878e-17, 0, 0, 0, 1, 0, 0, 0, 0, -0.00613319, -0.0113985, 0, 0,
      0.0211324, 0, 0, 0, -0.0700703, 0, 1.249e-16, 0, 0, 0, -0.0210687, 0, 0,
      0, 1.11022e-16, 0, 0, 0, 1, 0.00194066, 0, 0, 0, -0.0213299, 0, 0, 0,
      -0.0185897, 0, 0, 0, 0.0703197, 0.0237511, 0, 0, 0, -0.0176617, 0, 0, 0,
      0.0148613, 0, 0, 0, 0.00194066, 1, 0, 0, 0, 0.00980207, 0, 0, 0,
      4.12241e-16, 0, 0, 0, 1.249e-16, 0, 0.00415193, 0, 0, 0, -0.00971203, 0,
      0, 0, 0.00833903, 0, 0, 0, 0, 1, 0, 0, 0, 2.88549e-16, 0, 0, 0,
      3.43734e-16, 0, 0, 0, 0, 0, -0.0433016, 0, 0, 0, 0.0310482, 0, 0, 0,
      -0.0111976, 0, 0, 0, 0, 1, 0, 0, 0, 9.21638e-18, 0, 0, 0, 4.41139e-16, 0,
      0, 0, 0, 0, -0.0138466, -0.00256163, 0, 0, 0.0185842, 0, 0, 0,
      -0.00613319, 0, 0, 0, 0, 1, 0.00473485, 0, 0, -8.01766e-17, 0, 0, 0,
      -2.61143e-16, 0, -0.0176617, 0, 0, -0.00256163, 0.05203, 0, 0, 0.00561224,
      -0.00671312, 0, 0, -0.0113985, -0.0213299, 0.00980207, 0, 0, 0.00473485,
      1.19531, 0, 0, -0.0103735, 0.0155603, 0, 0, 0.0210687, -0.0210687, 0,
      -0.00971203, 0, 0, 0, 0.0246417, 0, 0, 0, -0.0263528, 0, 0, 0, 0,
      2.88549e-16, 0, 0, 0, 1, 0, 0, 0, 1.2767e-16, 0, 0, 0, 0, 0, 0.0310482, 0,
      0, 0, -0.0296094, 0, 0, 0, 0.0227903, 0, 0, 0, 0, 9.21638e-18, 0, 0, 0, 1,
      0, 0, 0, 4.28107e-17, 0, 0, 0, 0, 0, 0.0185842, 0.00561224, 0, 0,
      -0.030451, 0, 0, 0, 0.0211324, 0, 0, 0, 0, -8.01766e-17, -0.0103735, 0, 0,
      1, 0, 0, 0, -6.61878e-17, 0, 0.0148613, 0, 0, 0, -0.00671312, 0, 0, 0,
      0.0281051, 0, 0, 0, -0.0185897, 4.12241e-16, 0, 0, 0, 0.0155603, 0, 0, 0,
      1, 0, 0, 0, 1.11022e-16, 0, 0.00833903, 0, 0, 0, -0.0263528, 0, 0, 0,
      0.0451492, 0, 0, 0, 0, 3.43734e-16, 0, 0, 0, 1.2767e-16, 0, 0, 0, 1, 0, 0,
      0, 0, 0, -0.0111976, 0, 0, 0, 0.0227903, 0, 0, 0, -0.0414569, 0, 0, 0, 0,
      4.41139e-16, 0, 0, 0, 4.28107e-17, 0, 0, 0, 1, 0, 0, 0, 0, 0, -0.00613319,
      -0.0113985, 0, 0, 0.0211324, 0, 0, 0, -0.0700703, 0, 0, 0, 0,
      -2.61143e-16, 0.0210687, 0, 0, -6.61878e-17, 0, 0, 0, 1, 0, 0.00194066, 0,
      0, 0, -0.0213299, 0, 0, 0, -0.0185897, 0, 0, 0, 0.0703197, 1.249e-16, 0,
      0, 0, -0.0210687, 0, 0, 0, 1.11022e-16, 0, 0, 0, 1;

  bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_overlap, 1);
  if (!check_overlap) {
    cout << "ref" << endl;
    cout << overlap_ref << endl;
    cout << "result" << endl;
    cout << overlap.Matrix() << endl;
  }

  AOCoulomb coulomb;
  coulomb.Fill(auxbasis);
  Eigen::MatrixXd coulomb_ref =
      Eigen::MatrixXd::Zero(auxbasissize, auxbasissize);
  coulomb_ref << 0.177458, 0, 0, 0, 0.00173945, 0, 0, 0, -5.0899e-16, 0, 0, 0,
      -6.59195e-16, 0.00737579, 0, 0, 0, -0.0046258, 0, 0, 0, 0.00251989, 0, 0,
      0, 0.00101389, 0, 0.177458, 0, 0, 0, 5.72107e-16, 0, 0, 0, 2.60091e-16, 0,
      0, 0, 0, 0.00186685, 0, 0, 0, -0.00372032, 0, 0, 0, 0.00167281, 0, 0, 0,
      0, 0, 0.177458, 0, 0, 0, -4.81256e-16, 0, 0, 0, 4.5071e-16, 0, 0, 0, 0,
      -0.0128239, 0, 0, 0, 0.00741287, 0, 0, 0, -0.000167197, 0, 0, 0, 0, 0,
      0.177458, 0.000840236, 0, 0, 1.46075e-16, 0, 0, 0, 5.76574e-17, 0, 0, 0,
      0, -0.00579144, -0.000946348, 0, 0, 0.00629574, 0, 0, 0, -9.15777e-05, 0,
      0.00173945, 0, 0, 0.000840236, 0.659262, 0, 0, -0.00184086, 0.0027613, 0,
      0, 0.00373881, -0.00373881, -0.0046258, 0, 0, -0.000946348, 0.267231, 0,
      0, 0.00207334, -0.000303246, 0, 0, -0.00421098, -0.00776951, 0,
      5.72107e-16, 0, 0, 0, 0.177458, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0,
      -0.00372032, 0, 0, 0, 0.00868561, 0, 0, 0, -0.00648608, 0, 0, 0, 0, 0,
      -4.81256e-16, 0, 0, 0, 0.177458, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0,
      0.00741287, 0, 0, 0, -0.0078581, 0, 0, 0, 0.00504498, 0, 0, 0, 0, 0,
      1.46075e-16, -0.00184086, 0, 0, 0.177458, 0, 0, 0, -1.11731e-17, 0, 0, 0,
      0, 0.00629574, 0.00207334, 0, 0, -0.0104345, 0, 0, 0, 0.00329104, 0,
      -5.0899e-16, 0, 0, 0, 0.0027613, 0, 0, 0, 0.177458, 0, 0, 0, 1.74513e-15,
      0.00251989, 0, 0, 0, -0.000303246, 0, 0, 0, 0.00930186, 0, 0, 0,
      -0.00349634, 0, 2.60091e-16, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0.177458, 0,
      0, 0, 0, 0.00167281, 0, 0, 0, -0.00648608, 0, 0, 0, 0.0142938, 0, 0, 0, 0,
      0, 4.5071e-16, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0.177458, 0, 0, 0, 0,
      -0.000167197, 0, 0, 0, 0.00504498, 0, 0, 0, -0.0138951, 0, 0, 0, 0, 0,
      5.76574e-17, 0.00373881, 0, 0, -1.11731e-17, 0, 0, 0, 0.177458, 0, 0, 0,
      0, -9.15777e-05, -0.00421098, 0, 0, 0.00329104, 0, 0, 0, -0.0183512, 0,
      -6.59195e-16, 0, 0, 0, -0.00373881, 0, 0, 0, 1.74513e-15, 0, 0, 0,
      0.177458, 0.00101389, 0, 0, 0, -0.00776951, 0, 0, 0, -0.00349634, 0, 0, 0,
      0.0190278, 0.00737579, 0, 0, 0, -0.0046258, 0, 0, 0, 0.00251989, 0, 0, 0,
      0.00101389, 0.177458, 0, 0, 0, 0.00173945, 0, 0, 0, -5.0899e-16, 0, 0, 0,
      -6.59195e-16, 0, 0.00186685, 0, 0, 0, -0.00372032, 0, 0, 0, 0.00167281, 0,
      0, 0, 0, 0.177458, 0, 0, 0, 5.72107e-16, 0, 0, 0, 2.60091e-16, 0, 0, 0, 0,
      0, -0.0128239, 0, 0, 0, 0.00741287, 0, 0, 0, -0.000167197, 0, 0, 0, 0,
      0.177458, 0, 0, 0, -4.81256e-16, 0, 0, 0, 4.5071e-16, 0, 0, 0, 0, 0,
      -0.00579144, -0.000946348, 0, 0, 0.00629574, 0, 0, 0, -9.15777e-05, 0, 0,
      0, 0, 0.177458, 0.000840236, 0, 0, 1.46075e-16, 0, 0, 0, 5.76574e-17, 0,
      -0.0046258, 0, 0, -0.000946348, 0.267231, 0, 0, 0.00207334, -0.000303246,
      0, 0, -0.00421098, -0.00776951, 0.00173945, 0, 0, 0.000840236, 0.659262,
      0, 0, -0.00184086, 0.0027613, 0, 0, 0.00373881, -0.00373881, 0,
      -0.00372032, 0, 0, 0, 0.00868561, 0, 0, 0, -0.00648608, 0, 0, 0, 0,
      5.72107e-16, 0, 0, 0, 0.177458, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0, 0,
      0.00741287, 0, 0, 0, -0.0078581, 0, 0, 0, 0.00504498, 0, 0, 0, 0,
      -4.81256e-16, 0, 0, 0, 0.177458, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0, 0,
      0.00629574, 0.00207334, 0, 0, -0.0104345, 0, 0, 0, 0.00329104, 0, 0, 0, 0,
      1.46075e-16, -0.00184086, 0, 0, 0.177458, 0, 0, 0, -1.11731e-17, 0,
      0.00251989, 0, 0, 0, -0.000303246, 0, 0, 0, 0.00930186, 0, 0, 0,
      -0.00349634, -5.0899e-16, 0, 0, 0, 0.0027613, 0, 0, 0, 0.177458, 0, 0, 0,
      1.74513e-15, 0, 0.00167281, 0, 0, 0, -0.00648608, 0, 0, 0, 0.0142938, 0,
      0, 0, 0, 2.60091e-16, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0.177458, 0, 0, 0, 0,
      0, -0.000167197, 0, 0, 0, 0.00504498, 0, 0, 0, -0.0138951, 0, 0, 0, 0,
      4.5071e-16, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0.177458, 0, 0, 0, 0, 0,
      -9.15777e-05, -0.00421098, 0, 0, 0.00329104, 0, 0, 0, -0.0183512, 0, 0, 0,
      0, 5.76574e-17, 0.00373881, 0, 0, -1.11731e-17, 0, 0, 0, 0.177458, 0,
      0.00101389, 0, 0, 0, -0.00776951, 0, 0, 0, -0.00349634, 0, 0, 0,
      0.0190278, -6.59195e-16, 0, 0, 0, -0.00373881, 0, 0, 0, 1.74513e-15, 0, 0,
      0, 0.177458;
  bool check_coulomb = coulomb.Matrix().isApprox(coulomb_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_coulomb, 1);
  if (!check_coulomb) {
    cout << "ref" << endl;
    cout << coulomb_ref << endl;
    cout << "result" << endl;
    cout << coulomb.Matrix() << endl;
  }

  AOKinetic kinetic;
  kinetic.Fill(dftbasis);

  Eigen::MatrixXd kinetic_ref =
      Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);

  kinetic_ref << 29.9595, 0, 0, 0, 1.03944e-15, 0, 0, 0, -1.98878e-16,
      -0.011401, 0, 0, 0, -0.000410714, 0, 0, 0, 0.0614917, 0, 29.9595, 0, 0, 0,
      1.19856e-15, 0, 0, 0, 0, -0.0116766, 0, 0, 0, 0.0272569, 0, 0, 0, 0, 0,
      29.9595, 0, 0, 0, -1.87629e-15, 0, 0, 0, 0, 0.0573089, 0, 0, 0,
      -0.0337547, 0, 0, 0, 0, 0, 29.9595, 0, 0, 0, 3.30788e-16, 0, 0, 0, 0,
      0.0275401, 0, 0, 0, -0.0225031, 0, 1.03944e-15, 0, 0, 0, 29.9595, 0, 0, 0,
      -7.11277e-16, -0.000410714, 0, 0, 0, 0.012974, 0, 0, 0, -0.0734951, 0,
      1.19856e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0, 0.0272569, 0, 0, 0, -0.0734894,
      0, 0, 0, 0, 0, -1.87629e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0, -0.0337547, 0,
      0, 0, 0.0488035, 0, 0, 0, 0, 0, 3.30788e-16, 0, 0, 0, 29.9595, 0, 0, 0, 0,
      -0.0225031, 0, 0, 0, 0.0785723, 0, -1.98878e-16, 0, 0, 0, -7.11277e-16, 0,
      0, 0, 29.9595, 0.0614917, 0, 0, 0, -0.0734951, 0, 0, 0, 0.0237354,
      -0.011401, 0, 0, 0, -0.000410714, 0, 0, 0, 0.0614917, 29.9595, 0, 0, 0,
      1.03944e-15, 0, 0, 0, -1.98878e-16, 0, -0.0116766, 0, 0, 0, 0.0272569, 0,
      0, 0, 0, 29.9595, 0, 0, 0, 1.19856e-15, 0, 0, 0, 0, 0, 0.0573089, 0, 0, 0,
      -0.0337547, 0, 0, 0, 0, 29.9595, 0, 0, 0, -1.87629e-15, 0, 0, 0, 0, 0,
      0.0275401, 0, 0, 0, -0.0225031, 0, 0, 0, 0, 29.9595, 0, 0, 0, 3.30788e-16,
      0, -0.000410714, 0, 0, 0, 0.012974, 0, 0, 0, -0.0734951, 1.03944e-15, 0,
      0, 0, 29.9595, 0, 0, 0, -7.11277e-16, 0, 0.0272569, 0, 0, 0, -0.0734894,
      0, 0, 0, 0, 1.19856e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0, 0, -0.0337547, 0,
      0, 0, 0.0488035, 0, 0, 0, 0, -1.87629e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0,
      0, -0.0225031, 0, 0, 0, 0.0785723, 0, 0, 0, 0, 3.30788e-16, 0, 0, 0,
      29.9595, 0, 0.0614917, 0, 0, 0, -0.0734951, 0, 0, 0, 0.0237354,
      -1.98878e-16, 0, 0, 0, -7.11277e-16, 0, 0, 0, 29.9595;
  bool check_kinetic = kinetic.Matrix().isApprox(kinetic_ref, 0.00001);

  BOOST_CHECK_EQUAL(check_kinetic, 1);
  if (!check_kinetic) {
    cout << "ref" << endl;
    cout << kinetic_ref << endl;
    cout << "result" << endl;
    cout << kinetic.Matrix() << endl;
  }
}

BOOST_AUTO_TEST_SUITE_END()
